clear; clc;
% 扫描 AE, 在第三方程加上 AE*cos(0.6*t) 的周期电场

% 参数设定
a     = 0.7; 
b     = 0.8; 
c     = 0.1; 
alpha = 0.25; 
beta  = 0.01; 
delta = 0.1;

% 初始状态
y0 = [0.2; 0; 0.01];

% 积分时间
T_total = 1000;   % 总模拟时间
T_trans = 200;    % 舍去瞬态时间
dt = 0.01;       % 用于 LE 估计的重新归一化间隔
h = 0.01;        % RK4 步长

% 扫描 AE 范围
AE_range = linspace(0, 4.5, 700);

% 预分配
bifurAE = [];     % 分岔图数据 [AE, x_peak]
LE_vals = zeros(size(AE_range));

% 主要驱动 us=0.37*cos(0.7*t) 的周期
T_period = 2*pi / 0.7;

for i = 1:length(AE_range)
    AE_val = AE_range(i);

    params = struct('a', a, 'b', b, 'c', c, 'alpha', alpha, ...
                    'beta', beta, 'delta', delta, 'AE', AE_val);

    % 数值积分：调用 RK4 与含 E 的模型 hybridNeuron_AE
    [t, Y] = rk4(@(tt,yy) hybridNeuron_AE(tt, yy, params), [0, T_total], y0, h);

    % 舍去瞬态数据
    idx = (t > T_trans);
    t_post = t(idx);
    x_post = Y(idx,1);

    % 使用 findpeaks 检测 x(t) 的局部极大值
    [pks, ~] = findpeaks(x_post, t_post);

    % 分岔图数据
    bifurAE = [bifurAE; repmat(AE_val, length(pks), 1), pks]; % 将所有峰值 pks 存入

    % 计算 LE
    LE_vals(i) = LE_AE(params, y0, T_total, dt, h);
end

% (a) 分岔图
subplot(1,2,1);
plot(bifurAE(:,1), bifurAE(:,2), 'r.', 'MarkerSize', 5);
xlabel('A_E');
ylabel('x_{peak}');
title('Bifurcation diagram vs A_E');

% (b) LE vs A_E
subplot(1,2,2);
plot(AE_range, LE_vals, 'g-', 'LineWidth', 1.5); hold on;
yline(0, 'r--', 'LineWidth', 1.2);
xlabel('A_E');
ylabel('LE');
title('LE vs A_E');